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Modern split component evaluations of the flux vector jacobians are thoroughly 
analyzed for equilibrium-gas average-state determinations. It is shown that all such 
derivations satisfy a fundamental eigenvalue consistency theorem. A conservation - 
variable average state is then developed for arbitrary equilibrium-gas equations of 
state and curvilinear-coordinate fluxes. 

Original expressions for eigenvalues, sound speed, Mach number, and eigenvec- 
tors are then determined for a general average jacobian, and it is shown that the 
average eigenvalues, Mach number, and eigenvectors may not coincide with their 
classical pointwise counterparts. 

A general equilibrium-gas equation of state is then discussed for conservation- 
variable CFD Euler formulations. The associated derivations lead to unique com- 
patibility relations that constrain the pressure jacobian derivatives. Thereafter, 
alternative forms for the pressure variation and average sound speed are developed 
in terms of two average pressure jacobian derivatives. Significantly, no additional 
degree of freedom exists in the determination of these two average partial deriva- 
tives of pressure. Therefore, they are simultaneously computed exactly without 
any auxiliary relation, hence without any geometric solution projection or arbitrary 
scale factors. 

Several alternative formulations are then compared and key differences high- 
lighted with emphasis on the determination of the pressure variation and average 
sound speed. The relevant underlying assumptions axe identified, including some 
subtle approximations that are inherently employed in published average-state pro- 
cedures. 

Finally, a representative test case is discussed for which an intrinsically exact 
average state is determined. This exact state is then compared with the predictions 
of recently published methods, and their inherent approximations are appropriately 
quantified. 
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Introduction 

Since the notable paper by Roe [3], the classical average state and associated 
flux difference splitting algorithm have attracted the attention of several CFD re- 
searchers. Liou and van Leer [5] use this approach, and compare it to several flux 
vector splitting schemes for selected quasi 1-D problems, while Simpson [6] conducts 
a similar analysis for various 2-D airfoil flow field determinations. One often cited 
advantage of an average-state flux difference splitting algorithm lies in its reduced 
amount of intrinsic dissipation with respect to that inherent in the classical Steger- 
Warming and van Leer flux vector splittings. Furthermore, its capability to yield 
crisp shock resolutions makes it an attractive algorithm. These attributes have thus 
promoted continued research on real gas extensions, which ushered in the conspic- 
uous works of Glaister [4], Liou, van Leer &; Shuen [8], and more recently Vinokur 
& Montagne [9], among others. 

The linking thread of these fundamental contributions is the utilization of the 
classical ideal-gas average state as a building block in seeking a kindred real gas 
average state. As detailed in Section 4, all such derivations separate the pressure 
from the kinetic/energy terms, which corresponds to a decomposition of the flux 
vector into convection and acoustic- wave components. The kinematic/energy rela- 
tions are then exactly satisfied by the ideal-gas average state, whereas the pressure 
variation relation is evaluated at some specially devised average partial derivatives 
of pressure. The advantage of this procedure is that for an n-dimensional formula- 
tion, one eigenvalue of the associated average Jacobian coincides with the ideal-gas 
average velocity, with algebraic and geometric multiplicity equal to n, while the 
other eigenvalues resemble the eigenvalues of the pointwise Jacobian. 

In the published procedures [4,8,9], the terminal average Jacobian is obtained 
by restating the pressure Jacobian in terms of two average thermodynamic pressure 
derivatives. The pivotal equation to be satisfied is then the pressure variation 
relation in terms of these two average derivatives, which are essentially regarded as 
two auxiliary independent variables. Hence, the pressure variation relation is viewed 
as a single admissible-state straight line equation in the two initially unknown partial 
derivatives of pressure. Therefore, an absence on mathematical closure is perceived, 
and additional constraints are then devised to complement the pressure variation 
equation. As detailed in Section 9, these constraints are inherently based upon 
an ideal-gas average-state evaluation of selected expressions within the jacobian 
pressure variation, a choice of convenient integration paths and associated auxiliary 
states, and a utilization of geometric trial solution projections and scale factors. 

The specific avenue to develop the terminal average pressure derivatives then 
distinguishes each procedure. For example, in conjunction with an extra average 
mass-specific internal energy, Glaister [4] introduces additional intermediate states, 
whereas Liou et al. [8] utilize additional auxiliary partial derivatives of pressure. 
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These procedures actually correspond to an approximate integration of the pressure 
differential, upon selection of convenient integration paths, as discussed in Section 
9. This notion is expressly utilized by Vinokur et al. [9] who provide an exact 
integral representation for these derivatives. They then integrate these expressions 
numerically for an alternative determination of another set of average partial deriva- 
tives. These constructions employ an independent internal energy variation, even 
though this variation directly depends upon the state variable, as detailed in Sec- 
tion 8. Futhermore, this approach generally yields approximate solutions that may 
not satisfy the pressure variation relation. Consequently, the terminal derivatives 
that do satisfy this relation are determined by projecting the approximate solutions 
onto the admissible-state straight line, according to a set of scale factors. 

Putting aside the legitimate concern that additional intermediate states may 
not be permissible by the given real gas equations of state, the crucial issue is 
that these procedures require several assumptions and introduce different sets of 
average partial derivatives of pressure and sound speed expressions. And even 
though these relations satisfy a fundamental eigenvalue consistency theorem, they 
generally do not coincide with one another, or they may not exist in the domain of 
the analytical partial derivatives, or they may correspond to conservation-variable 
evaluations states outside the hyper-segment connecting the two given cell states, 
or they may yield indeterminate expressions. Thus, they may correspond to a 
state weakly related to the given cell states, notwithstanding the exact algebraic 
enforcement of the pressure variation expression. 

In this connection, a representative test case in Section 10 illustrates the differ- 
ence between an exact average state and the predictions of some of these methods. 
The average-state jacobian may thus become less representative of the local wave 
interaction patterns, thereby defeating the concept of an accurate approximate local 
Riemann solver for which the average-state notion was introduced in the first place. 
Moreover, citing Vinokur [9], the amount of intrinsic dissipation in these conditions 
may become a vaguely identified quantity. 

In this report it is shown that the published equilibrium-gas average states [4,8,9] 
satisfy a fundamental eigenvalue consistency theorem. According to this theorem, 
if the flow is locally subsonic, the eigenvalue set of any associated average jacobian 
will contain mixed sign eigenvalues, whereas for a locally supersonic flow all average 
eigenvalue sets will contain uniform sign eigenvalues. The most recent average- 
state determination procedures are discussed, and a set of conservation- variable 
average states is then developed for arbitrary equilibrium-gas equations of state 
and curvilinear- coordinate fluxes. These states correspond to an internal average 
state belonging to the hyper-segment connecting the given computational cell left 
and right states. 

A detailed characteristic analysis of the general curvilinear coordinate average 
jacobian is performed which leads to original expressions for eigenvalues, sound 
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speed and Mach number in terms of the pressure jacobian. These relations, which 
revert to the classical thermodynamic expressions, show that the average eigenval- 
ues, Mach number and eigenvectors may not coincide with the classical pointwise 
forms. 

Since for an equilibrium-gas formulation, pressure is a function of only two ther- 
modynamic variables, it is then shown that the pressure jacobian derivatives are not 
functionally independent. On the contrary, they satisfy a compatibility constraint 
that leads to alternative expressions for the pressure variation and average sound 
speed in terms of only two average pressure jacobian derivatives. These expres- 
sions do not explicitly depend upon the thermodynamic internal energy, since in a 
CFD formulation this variable is not independent, but is expressed in terms of the 
conservation state variable. Considering that within the mean value theorem reso- 
lution of the pressure variation, the pressure derivatives are intrinsic functions of a 
conservation state variable, it is shown that no additional degree of freedom exists 
in their computation. Therefore, these derivatives are simultaneously determined 
exactly without resorting to any geometric solution projection and arbitrary scale 
factors. The developed formulae, which bear a striking similarity with the classi- 
cal thermodynamic relations, apply for totally arbitrary equilibrium-gas equations 
of state, and yield as approximate cases several expressions utilized in reported 
procedures. These alternative formulations are then compared and their intrinsic 
underlying assumptions are identified including some inherently employed subtle 
approximations. 

To conclude, a representative test case is analyzed for which an intrinsically 
exact average state is determined without separating the kinematic/energy terms 
from the pressure relation in the flux difference resolution. Using this average state, 
the related exact average thermodynamic pressure derivatives and sound speed are 
computed. These results are then compared with the predictions of recent methods 
and their inherent approximations are appropriately quantified. 

This report is organized in 10 Sections. Section 1 delineates the relevant form 
of the curvilinear- coordinate Euler equations, while Section 2 contains several the- 
oretical average state considerations. Section 3 presents the eigenvalue consistency 
theorem, while Section 4 details the split component evaluation of the flux vector 
jacobian. Sections 5-6 axe devoted to the development of the conservation-variable 
average state for one- and multi- dimensional curvilinear coordinate fluxes, and a 
discussion of the ideal gas average state. Section 7 details the characteristic analysis 
derivation of the eigenvalues and the sound speed expression in terms of the pres- 
sure jacobian, while the pivotal average partial derivatives of pressure are directly 
determined in Section 8. Finally, several comparisons with other average-state for- 
mulations axe presented in Section 9, while representative results are discussed in 
Section 10. 
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1 Governing Equations 


The cartesian coordinate Euler equations in conservation law form are 


9q dfj 
dt dx; 


= 0 


on ^ + xfi, t > t 0 , t,t 0 e$l + , ft C £ n , 1 < n < 3 (1) 

where n denotes the spatial dimensionality. The state variable q = q(x,t), and the 
flux vector fj = f } {q) are 
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where p is the fluid density, m = pu is the linear momentum vector, p is the static 
pressure, 6j = {£!•}, 1 < i < n, is the Kronecker delta mixed tensor, and E is the 
volume specific total energy. System (l)-(2) is then closed with an equation of state 
for pressure. For a homogeneous gas in thermo-chemical equilibrium, this equation 
of state may be expressed as a function of q, differentiable almost everywhere [1] in 
ft, of the form 

p = p(q(x,t)) (3) 

as detailed in Section 8. 

With these specifications, the conservation state variable q constitutes the differ- 
etiation variable for determining the jacobian matrix Aj(q) of fj(q)- Consequently, 
any other variable introduced to compute Aj(q) will depend upon q. Hence, the 
partial derivatives of p(q) in Aj(q) coincide with the jacobian of p(q) with respect 
to q. This pressure jacobian can certainly be computed by way of the differenti- 
ation chain rule and the thermodynamic derivatives of pressure. Nonetheless, the 
required additional thermodynamic variables axe themselves functions of q. These 
considerations have significant repercussions on the interpretation and evaluation 
of the pressure variation relation, as thoroughly documented in Section 8. 

For arbitrary-geometry computations, the cartesian coordinates x are trans- 
formed into curvilinear coordinates rj via the relation 

* = x(rj ) (4) 

Function (4) must be single valued and have continuous first partial derivatives. 
The jacobian J of the coordinate transformation (4) is 


The relation between the Xj and r] 3 partial derivatives is then 


dr)k 


dr) k \\dxjf \dxjj 


(6) 


Therefore, the partial derivatives with respect to the cartesian coordinates x 3 are 
expressed using the inverse of (6) as 


d'Xj 


= J 


-l 


drj k 


= [det J -J- 1 ] {^-1 

det J J \ dr)k J 


(7) 


where det J does not vanish and keeps a uniform sign in ft. With these specifica- 
tions, the matrix [det J • J -1 ] is denoted as \e 3k ). The metric data e jk , with con- 
tinuous partial derivatives in their definition domain, then satisfy the fundamental 
invariance relation 

E^=0 , \<j<n (8) 

tx dr )k 

Inserting (7)-(8) into (1) yields the contravariant conservation law system 


£? , 1 9fk _ 

dt det J dr]k 

where the contravariant flux /* is defined as 


(9) 


e jk mj 
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fk = e jk {//} = e jk fj = e jk -£m + e }kP 6 3 

e lk-^{E + p) 


( 10 ) 


The curvilinear coordinate average states are then determined for the Jacobian of 
this contravariant flux, as detailed in Section 4. 


2 Theoretical Average States 


A set of average states for expressing the flux variation exists by virtue of the mean 
value theorem. For an ideal gas, the original resolution [3] accrues from one average 
state and one average enthalpy. For a real gas, an analogous resolution also depends 
upon the partial derivatives of pressure. 

The variation of each flux component fj in (2), (10) can be expressed by way of 
the multi-dimensioned mean value theorem as 





(,* +a2 


(ii) 
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where the partial derivatives in ( 11 ) are evaluated at the internal state 



( 12 ) 


In ( 11 )-( 12 ) superscripts L and R denote the available computational-cell left and 
right states respectively. Furthermore, each weight parameter a* satisfies the rela- 
tion 

0.0 < a e < 1.0 , l<£<n + 2 (13) 

In general, these parameters may not coincide with one another, and hence the 
internal states < 7 ^, 1 < £ < n + 2 , may differ from one another. Therefore, the 
variation of the flux vector ( 2 ), ( 10 ) is cast as 

A] (<?\ 9 2 , • ■ • , <T +2 ) Aq = Afj (14) 

where, following (11), the i th row of Aj is evaluated at q* . Furthermore, more than 
one state may exist for a fixed i . Thus, the average state is in general non-unique. 

While the mean value theorem evaluation relies on at most n + 2 states, the most 
recent developments only need fewer such states upon introduction of an additional 
average enthalpy. This reduction accrues from an independent application of the 
mean value theorem to specific components in the resolution of each flux vector 
component variation. For example, (11) can be cast as 


A/? = A ff + A ff = 


dq 


(? L + * {q R - /)) Aq + (« L + * (g* - , L )) A , 

( 15 ) 

whereupon, the mean value theorem independently applied to A /J and A /“ yields 


A ff = (q L + S‘- (q R ~ q L )) Aq (16) 

A ff = < 17 > 

As shown in Section 4, this procedure corresponds to a convection/acoustic- wave 
flux vector decomposition that is intrinsically utilized in [4,8,9]. 

Even though several CFD papers [3, 5, 8 , 9] employ the appealing notation 


Aj (?) Aq = A fj 


(IS) 


which relies on a free interpretation of the multi-dimensional mean value theorem 
and presumes the existence of a single average state q valid for each and every flux 
vector component, the original ideal-gas average state [3] violates the mean value 
theorem (11) in the flux vector energy component variation, as detailed in Section 6 . 
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Nevertheless, the original average state successfully yields a closed form resolution 
for (11) using an additional average enthalpy H . Thus (18) can be more accurately 
cast as 

Aj (q,H) Aq = (19) 

and is at variance with the developments in [3]. As proved in Section 5, expression 
(19) for an ideal gas can be alternatively cast as 

Aj(q,p)&q = Afj (20) 


upon identification of an internal-average total energy with associated average pres- 
sure p. For a general equilibrium gas equation of state, expression (20) is expanded 


as 



m,E 



A q = A/, 


( 21 ) 


r,m/ 


where a tilde denotes average partial derivatives. In this form, relation (21) indi- 
cates the crucial functional dependence of the flux vector jacobian upon the partial 
derivatives of pressure. 


3 Eigenvalue Invariance Relations 

Considering the non-uniqueness of the average states, several dissimilar average ja- 
cobians can be developed for a single state- variable pair. Consequently, the spectra 
of these jacobians will not coincide with one another. While this situation may be 
unavoidable, physical consistency demands that, for a local flow of a specific charac- 
ter (subsonic or supersonic), all associated average jacobians for any average state 
must correspond to a uniform set of eigenvalues. Therefore, if the flow is locally 
subsonic, the eigenvalue set of any associated average jacobian must contain posi- 
tive and negative eigenvalues, whereas for a locally supersonic flow, the eigenvalue 
set of all associated average jacobians can only contain uniform sign eigenvalues. 
If this property were not always met, then for some subsonic flows, a flux differ- 
ence splitting implementation may lead to a fully upwind scheme for all flux vector 
components, which may in turn engender numerical solution instability. As proven 
by the following n-dimensional-formulation theorem, those average jacobians that 
satisfy (21) and share the ideal-gas average velocity eigenvalue, with algebraic and 
geomentric multiplicity equal to n, display consistent eigenvalue sets with mixed- 
or uniform-sign eigenvalues for a subsonic or supersonic flow state respectively. 

Definition : two non-singular matrices are defined eigenvalue consistent when 
their spectra have the same number of positive and negative eigenvalues. 
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Theorem 1 Let two non-singular and diagonalizable average jacobians correspond 
to any single set of state variable and flux variations A q and A/. These jacobians 
are then always eigenvalue consistent for supersonic flows. For subsonic flows they 
are eigenvalue consistent if for an n dimensional formulation they share one eigen- 
value with algebraic and geometric multiplicity equal to n. 

Proof . For any state variable variation A q, the resolution of the associ- 
ated flux variation A / for any two generally different average states can 
be expressed as 

A ( ?I ’ ft ’(ir),) Aq = Af ' 

Given the purely algebraic nature of the following derivations, ( 22 ) is 
abridged as 

A x Aq = A / , A 2 Aq = Af (23) 

Since each matrix in (23) is diagonalizable, each possesses a complete 
set of eigenvectors. Hence Aq can be expressed as 

A q = X\b , Aq = X 2 c (24) 

In (24), X\ and X 2 are the eigenvector matrices for A\ and A 2 respec- 
tively, while b and c are the associated eigenvector coefficient set. From 
(24) it follows 

X 1 b = X 2 c (25) 

Furthermore, expressing the matrices A x and A 2 by way of a similarity 
transformation yields 

X 1 A 1 Xf 1 Aq = Af , X 2 A 2 Xf 1 Aq = Af (26) 

where Ai and A 2 are the diagonal matrices containing the eigenvalues 
of A\ and A 2 respectively. Inserting (24) into (26) leads to 


Aq = Af (22) 


X 1 A 1 b = Af , X 2 A 2 c = Af (27) 

whence 

X 1 A 1 b = X 2 A 2 c (28) 

Thereafter, (25) and (28) lead to the following expressions for b 

b = (aT 1 - 1 X 2 )c (29) 

b = A 7 1 A 2 c (30) 
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It must be shown that when Aj displays a certain number of positive and 
negative eigenvalues, A 2 has the same number of positive and negative 
eigenvalues. Therefore, for any A 5, hence for any 6, consider the positive 
bilinear form 

0 < ^IiAi b (31) 

In (31), Xj is a diagonal matrix with entries a{, 1 < j < n + 2, equal 
to +1 or —1 only, depending on the algebraic sign of the eigenvalues in 
Ai, such that the matrix X1A1 is positive definite. Expressing b T and b 
in (31) by way of (29) and (30) respectively yields 

0 < c T (X^X 2 ) T l x (Af 1 ^) A 2 c (32) 

According to (31), the matrix in (32) is positive definite. Consequently, 
all of its principal minors 

det [M\i = det [(A 1 - 1 X 2 ) T Jj (Xf 1 X?) A 2 ] (33) 

must be positive [2]. Considering the diagonality of A 2 , this matrix 
minor positivity property is expressed as 

0 < det (X^X 2 ) T I 1 (X{ 1 X 2 ) A 2 ] 

= det [(Xr 1 ^) 7 !! (Xr x X 2 )l • nA J 2 , l < i < n + 2 (34) 

where A 2 , 1 < j < n + 2, denote the eigenvalues in A 2 . For a supersonic 
flow, all the eigenvalues in Aj have the same algebraic sign. Hence, all 
of the entries o\ in T\ are equal to one another, and (34) thus becomes 

0 < det [(X 1 - 1 X 2 ) T (Xf : l X 2 )| • J] H A 2) , 1 < * < n + 2 (35) 

>- J ‘ j = 1 

Furthermore, the matrix 

Y = (Xf 1 X 2 ) T (Xr x X 2 ) (36) 

is certainly positive definite, as well as symmetric, since it corresponds 
to the matrix product of a non-singular matrix and its transpose. Con- 
sequently, the matrix minors 

det [F],. = det (Xr 1 X 2 ) T (x i - 1 X 2 )] , l<i<n + 2 (37) 
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are also positive. Therefore, it follows from (35) that 

> 0 , 1 < j < n + 2 (38) 

Hence, if Aj > 0, 1 < j < n + 2, then a\ — +1, 1 < j < n + 2, and 
A^ > 0, 1 < j < n + 2, whereas if Aj < 0, 1 < j < n + 2, then crj = —1, 

1 < j < n + 2, and \[ < 0, 1 < j < n + 2, which proves eigenvalue 
consistency. 

For subsonic flows, Ai has mixed sign eigenvalues. It thus has two eigen- 
values of opposite sign in addition to one eigenvalue with multiplicity n. 
Therefore, X\ has entries with different algebraic sign and consequently 
its determinant can be either positive or negative. For i = n + 2, ex- 
pression (34) yields 

n+2 

0 < det [Y] • det pi] • JJ A j 2 (39) 

j=i 

Consequently, 

n+2 

0 < det pi] • n < 4 °) 

j=i 

The eigenvalues A 2 , 1 < j < n + 2, cannot all have the same algebraic 
sign, otherwise, according to the previous conclusion, the same would 
hold true for Aj, 1 < j < n + 2, and the flow would be supersonic. 
Therefore, A 2 will display at least two eigenvalues of opposite sign. By 
hypothesis, Aj and A 2 already share one eigenvalue with multiplicity n, 
and for one- and three- dimensional formulations, n is an odd number. 

In this case, if the common eigenvalue is positive, respectively negative, 
det pj] is negative, respectively positive. Thus, from (40), det [A 2 ] is 
negative, respectively positive, but in either case A 2 has n eigenvalues 
with the same sign and two eigenvalues with opposite sign. For a two- 
dimensional formulation, n is an even number and det pi] is always 
negative. Therefore, from (40) det [A 2 ] will always be negative, and 
thus A 2 has again n eigenvalues with the same sign and two eigenvalues 
with opposite sign. Hence, also for subsonic flows Aj and A 2 have the 
same number of positive and negative eigenvalues. 

Finally, since this proof is valid for any two average jacobians that satisfy 
the stated hypotheses, it remains valid for all such jacobians □. 

By virtue of this theorem, a flux- difference splitting procedure will always yield 
the appropriate upwind scheme, for both subsonic and supersonic flows, for any 
average state that satisfies (21), provided that the associated average jacobian has 
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one eigenvalue with algebraic and geometric multiplicity equal to n. Therefore, it 
is not strictly indispensable to upwind scheme appropriateness that the eigenvalues 
of all average jacobians, associated with a single local flow character, coincide with 
one another. However, the closer the average eigenvalue sets are in an appropriate 
norm to the pointwise set, the more accurate the average-state flux vector difference 
resolution becomes. The following theorem shows that this eigenvalue convergence 
property directly depends upon the eigenvector matrices, as can be easily imagined. 

Theorem 2 If the matrices A\ and A 2 in Theorem 1 share the same eigenvector 
matrix X , then they also share the same eigenvalues. 

Proof . Theorem 1 yields the relations 

b - (X^X 2 ) c , b — A, -1 (Xr x X 2 ) A 2 c (41) 

Since X 1= X 2 = X , (41) yields 

c = AJ' 1 A 2 c (42) 

Therefore, A\ and A 2 have the same eigenvalues □. 

Consequently, it becomes convenient for the average eigenvector matrix to remain 
close in an appropriate norm to the pointwise eigenvector matrix. This is especially 
important for transonic or stagnation-point flows where one or more average eigen- 
values may vanish, hence some hypotheses of the eigenvalue consistency theorem 
may not be satisfied. 


4 Reference Jacobian and 

Current Evaluation Procedures 


Specific average states are determined such that the relations corresponding to the 
flux difference resolution 



(43) 


are formally satisfied. For a reference one- dimensional formulation, the average 
jacobian in (43) is 


m 2 / dp ' 

+ Wp, 


P 


m - _ m ( dp' 

P 2 P \dpj 


m } E 



Therefore, (43) yields the identity 



m R — m L — m R — m L 


(45) 


for the continuity equation, 



for the momentum equation, and 



for the energy equation. 

These three relations would constitute a very appealing formulation if a single 
average state q were efficiently identifiable for the evaluation of the jacobian state 
vector terms as well as pressure and its partial derivatives. Certainly (46)-(47) 
form a non-linear system of two equations for the three unknown components of 
q. Hence, the average state q could coincide with the solution of this system. 
However, no complete information exists to date regarding solution existence, let 
alone uniqueness. 

Therefore, the current formulations seek an average state by separating the 
pressure from the kinematic/energy terms in the momentum and energy equations, 
following (15)-(17). This procedure corresponds to a decomposition of the flux 
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vector fj(q) into the physically distinct convection and acoustic- wave components 
/'(g) and /“(g) 

Mi) = /;(«) + /*(?) ( 48 ) 

with definitions 




771 j 


P 

m 

{ E + p 


m = 



(49) 


The jacobian Aj(q) of fj(q ) is then 


, , , _ e/i _ a n , 


(50) 


Consequently, the resolution of the flux vector fj(q) can be expressed as in (15)-(17). 
The associated convection and acoustic-wave average states may then be distinct 
from one another considering the difference between these two processes. 

Thereafter, the variations in the kinematic/energy terms in the convection flux 
are expressed in terms of the ideal gas average state, whereas the pressure varia- 
tion in the acoustic- wave flux is evaluated by way of some specially devised average 
partial derivatives of pressure. The developed average density and velocity (or mo- 
mentum) that satisfy the momentum equation are then inserted into the energy 
equation. This technique yields an average enthalpy (or energy with associated 
pressure as shown in Section 5) in closed form which however does not satisfy the 
mean value theorem for the energy equation, as proven in Section 6. Nonetheless, 
a distinctive advantage of this modus operandi is that one eigenvalue of the as- 
sociated average jacobian is always the ideal-gas average velocity, with algebraic 
and geometric multiplicity equal to n, which satisfies the hypotheses of the eigen- 
value consistency theorem. Furthermore, the other eigenvalues resemble those of 
the pointwise jacobians, as detailed in Section 7. 

For the convection flux resolution, the separate evaluation method yields the 
identity 

( 51 ) 


m 


R 


m L = m R — m L 


for the continuity equation, 


(p R - p L ) + 2~ (m R - m L ) — 

K ) 2 

K)’ 



V / p \ / 

P 

P 


(52) 


for the momentum equation, and 
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for the energy equation. 

The terminal form of (54) is obtained by virtue of the pressure variation expres- 
sion associated with the acoustic flux. From (49)-(50) and for an n dimensional 
formulation, this expression is 



(55) 


which thus becomes the crucial relation to be satisfied. As fully detailed in Section 
9, the various procedures published thus far [4,8,9] differ in the specific avenue to 
evaluate this pressure relation, even though they share the common feature of alge- 
braically constructing sets of average partial derivatives of pressure assuming several 
geometric solution projection rules and associated scale factors. Consequently, the 
corresponding solutions remain inextricably tied up with these assumptions. On 
the contrary, the developments in Section 8 dispense with these assumptions by 
exactly determining the intrinsic acoustic-wave average state that satisfies (55). 

5 Conservation Variable Average States 

It is often reported that an average density p can be chosen freely [8], and that 
there is no physical basis for introducing specific expressions for p [9]. However, the 
classical average density p, as well as the average momentum m, and volume specific 
total energy E can all be analytically determined in a rigorous manner invoking the 
mean value theorem for the flux difference resolution (44). The resulting expressions 


p = p L + a (p R - p L ) (56) 

m = m L + a (rn R — m L ') (57) 

E = E l + a (E R - E l ) (58) 

p = p L + a (p R - p L ) (59) 


employ a single weight parameter a, 0 < a < 1, and apply for both one-dimensional 
and curvilinear-coordinate multi-dimensional formulations. 
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5.1 One-Dimensional Average States 

Expression (52) constitutes a quadratic equation for 


(!) 


( 60 ) 


with internal-state solution 


m 

rh _ p L 
P 




,R 


P L + \IP 


( 61 ) 


Furthermore, the insertion of expressions (56)-(57) into (61) yields a single equation 
for a with solution 


u 

+ < 
+ 

V 

\A 

h-* 

(62) 

P- \Jp L P R 

(63) 


With this solution, (56) becomes 


which coincides with the classical form. 

Inserting the average density and momentum into the two leading terms in the 
lhs of (53) yields the equality 


( 64 ) 


' m B(„R R m R m L \ 

~p E \P + ~ m )- e \Jr--^) 

Hence, (53) becomes 

which constitutes the average-state resolution of the flux of E. A similar operation 
sequence on (54) yields 

(m \ m R R m L , ^ ( m R m L \ m / R L \ 

A (7 P ) = j + J ' p ) <66) 

which constitutes the average-state resolution of the flux of p. Whereupon, inserting 
(61) into (65) yields 


E = E 


R 


+ e r - 


+ \fp R +\fp L 

= E L + a(E R -E L ) 


( 67 ) 
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where a is defined by (62). Hence, the average total energy E is also determined 
with exactly the same weight as the average density and momentum. Finally, the 
average pressure p is analogously expressed as 

p = p L + a(p R - p L ) (68) 

owing to the identical structure of (65) and (66). 


5.2 Multi-Dimensional Average States 

The derivations of the average state for the jacobian of the multi- dimensional flux 
(10) with embedded metric data tjk are similar to the previous developments, and 
show that the average one-dimensional conservation variables also apply for multi- 
dimensional curvilinear coordinate formulations. The variation of the contravariant 
flux fk is expressed as 


Ak 



Ag — A {ejkfj) 


(69) 


In order to simplify the notation, the partial derivatives of pressure in (69) are 
denoted as 


Pp = 


'dp 




PE 


m,E 

Therefore, the average state jacobian in (69) becomes 


■(») 


(70) 


p,m 


A k = A k (q,p,p p ,Pm,PE,ejk) = = 


e-tk 


0 


TJlj 


jk .~y TTli &ikPp 

_ P 


m , ,• /- \ m, 

r(® + J>) + «>»-* 


Pp 


m, rrij • __ 

(■fk ~ "f 6 ikPmi “I - &jk ~Oit i &ikPE 

P _ P _ 

/ - \ m, __ m, 

(£ + p) p mt , e jk —(l +p E ) 

v 7 P P 


^tk_ ( fk . 

P 


( 71 ) 

where, 1 < (i,^) < n are the row and column indices, for the n-dimensional res- 
olution of momentum. In the sequel, the various formulae are exemplified and 
demonstrated for a two-dimensional formulation. In this case, the average jacobian 
Ak becomes a 4 x 4 matrix, where the middle two rows and columns respectively 
derive from the middle row and column in (71) for = 1,2. For constant metric 
data ejk in each and every computational cell, the contravariant flux variation in 
(69) becomes 

A{ejkfj) = e jk A(fj) (72) 
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Therefore for the convection and acoustic- wave flux vector decomposition, expres- 
sion (69) yields the identity 

ei k (mf - mf) + e 2k (mf - m^) = e ]k mf - e jk mf (73) 


for continuity equation, 



eik ( P R ~ P L ) (75) 


for the longitudinal-momentum equation, 



e 2 t{p R ~P L ) (77) 


for the transverse-momentum equation, and 



for the energy equation, where (79) implies (75), and (77). 
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Expressing the principal-momentum flux variation in homogeneous form yields 



and since the metric data are linearly independent, the two expressions between 
brackets will simultaneously vanish. Whereupon, the first expression reverts to 
(60), hence the average density coincides with (56), while the average principal 
momentum is 

rhi ~ rnf + ot (mf — mf) (81) 

where a is defined by (62). Thereafter, inserting the average density and principal 
momentum into the second expression yields 

77^2 [—m[p R + mfp L ) + p (mi {rn R — = p L m R m R — p R m\m\ (82) 

which defines the average transverse momentum as 

ra 2 = m\ + a (m R — ) (83) 


where a is again defined by (62). The homogeneous form of the transverse - mo- 
mentum flux variation is then 



After interchanging the momentum subscripts, this expression reverts to the previ- 
ous principal-momentum flux variation, and consequently it is identically satisfied 
by the developed average density and momentum components. The average partial 
derivatives in expressions (75), (77) are then exactly determined in Section 8. 

The insertion of the average density and momentum components into the two 
leading terms in the lhs of (78) yields 





E 



ejk 



(85) 
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Hence, (78) becomes 

^ V J * p / ~ ,k p R E jk p 1 ' E \ ,fe p R jk p L ) + Cjfc p \ ^ / 

' ' ' . . ( 86 ) 
which constitutes the average-state resolution for the variation of the contravariant 

flux of E. A similar operation sequence on (79) yields 

/ m, \ mf B m f t ~f m ? m f\ m, / □ M 

r7 p J = Cife 7« p " = p _ eiA > ) + 6jk T ~ p ) 

(87) 

which constitutes the average-state resolution for the variation of the contravariant 
flux of p. Whereupon, inserting the average density and momentum components 
into (86) and considering the linear independence of the metric data yields the 
average volume specific total energy 

E = E l + o (E R - E l ) 

exactly as before. The average pressure is then again expressed as 

P = p L + a ( p R - p L ) 

owing to the identical structure of (86) and (87). 


( 88 ) 


(89) 


6 Reference Average State 

The original convection average state [3] relies on an average velocity and enthalpy, 
and can be expressed as 


P = yjp^p* ( 90 ) 



No average pressure p or volume specific total energy E is explicitly defined in the 
original derivations. However, the developed conservation variable average state 
(56)-(59) identically reverts to (90)-(92). This is proven by inserting the average 
density, momenta, energy and pressure states (56)-(59) into the definition of each 
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velocity component and enthalpy which yields 



(93) 

(94) 


Therefore, the developed convection average states (56)-(59) are totally equivalent 
to (90)-(92). It is occasionally reported that the average pressure should be obtained 
using its own equation of state evaluated at (90)-(92). However, if this course of 
action is selected, then the corresponding average total energy cannot coincide with 
the internal average (58). In either case, only in the enthalpy function do pressure 
and total energy emerge in the conservation variable jacobian. Therefore, any set 
of p and E that does not alter (94), will yield the same convection average state. 

As stated in Sections 2-3, the average enthalpy does not satisfy the mean value 
theorem for the energy equation. This is proved observing that the legitimate mean 
value theorem variables are the jacobian differentiation variables which coincide 
with the entries in q. Hence the total enthalpy 

#(<?) = + p{p,m,E)) (95) 

is itself a function of the dependent variable q. Consequently, its mean value theorem 
average is 

H mvt = ± (E + p(p,fh ,£)) (96) 

where superscript MVT denotes such an average. However, the average enthalpy 
(94) may not coincide with this expression. In fact, for the ideal gas equation of 
state (144), + he difference between (96) and (95) is 

(97) 

(98) 

where the terminal minimal variation form accrues from the internal average total 
energy (58). 

As detailed in Section 7, the average sound speed squared is exactly expressed 
as 
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and depends upon the average enthalpy (94). This dependence, which results from 
the insertion of the average density and momentum into the energy flux resolution 
(78), constitutes the fundamental reason why (99) may supply an external average 
sound speed, for some left and right states. For example, using (144), the ideal gas 
specification of the average sound speed squared becomes 


„2 _ 


= (7 - 1) 


771 


777 


2p 2 


+ (7 - 1) [H - 


777 


777 


( 100 ) 


The insertion of the average states (91) and (92) into (100) then yields 


a(a - 1)(7 - 1 ) / m R m 


c 2 = 


m 


R 


n>R 


^R 


(c L ) 2 + (l-a) (c R ) (101) 


which, as also determined in [7], adds an extra term to an internal average for c 2 . 
Conversely, if H in (99) were replaced by (96), this extra term would not emerge. 
However, (99) is the only correct expression compatible with the flux difference 
resolutions (43), (69). This result should not be considered detrimental, since the 
fundamental flux difference relations (43), (69) are rigorously satisfied for arbitrary 
left and right states. On the contrary, (99), hence (101), may be used to determine 
whether some external average eigenvalues occur in the computed flow field, and, 
consequently, whether a local enhancement of the basic flux difference scheme is 
warranted. 


7 Average Characteristic Relations 

A completely general eigenvalue analysis performed on the multi- dimensional curvi- 
linear - coordinate jacobian (71) yields an original expression for the speed of sound 
which is shown to revert identically to the classical thermodynamic relation. The 
same analysis yields the average eigenvalues, Mach number, and eigenvectors as 
functions of both the average state (56)-(59) and the partial derivatives of pressure. 
Importantly, these results apply for any specific form of these derivatives. 


7.1 Average Eigenvalues 


The trace of matrix (71), which coincides with the sum of its eigenvalues, is 


tr(A) = (» + 2)e* y + ^(^“) 


P,TTli,E,i^3 


m ? ( dp 

+e >‘fUi 


( 102 ) 


pjn 


Considering that in general the acoustic-wave average state corresponding to the 
average partial derivatives of pressure may differ from q , the trailing two-term ex- 
pression in (102) may not vanish since it may not satisfy the compatibility relation 
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(127) derived in Section 8. Hence, expression (102) may not reduce to the usual re- 
sult tr = ( n + 2)ejk Consequently, the eigenvalues of (71) may not coincide 

with the classical expressions, which for a three-dimensional cartesian-coordinate 
reference frame are ^ 

x _ \ - 

^* 1,3 ^ 5 ^* 4,5 

P 

where c denotes the speed of sound. In fact, the analytical solution of the charac- 
teristic equation 


m k 


± c 


(103) 


det 


Ak 




(104) 


yields 


m 7 


= e jk -A 

P 


— £jk -s ' + Cjk^dj i ^ (^jk^ dj “I - e jh e jk a 


1/2 


(105) 

(106) 


where Ud = {udj } r , 1 < J < n, and a 2 are directly obtained as part of the eigenvalue 
solution itself as 





dp \ m f dp\ E + p / dp \ 

dm ) e ,E P + \ d P/m.E + P \ 9E ),,m 


(108) 


With these expressions the sum of eigenvalues (105)-(106) then coincides with (102). 


7.2 Speed of Sound 

Relation (108) delivers an original alternative expression for the square of the speed 
of sound that identically reverts to the classical thermodynamic expression 


c 2 = 



(109) 


as directly obtained using the equation of state, the first principle of thermody- 
namics, and the definition of speed of sound. This coincidence accrues from a 
differentiation chain rule restatement of the partial derivatives of pressure in (108) 
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using expressions (124)-(126), detailed in Section 8. Insertion of these expressions 
into the generic form of (108) leads to 



which coincides with (109). Consequently, a 2 is indeed an alternative relation for c 2 . 
Whereupon, insertion into (108) of the pressure derivative compatibility relations 
(127), derived in Section 8, yields the average conservation- variable speed of sound 
squared 


c 


2 



m,E 



pjn 




an) 


which depends on only two partial derivatives, hence it closely resembles the classical 
thermodynamic relation (109). Expression (111) does not depend on an average 
internal energy, hence is not affected by any assumption connected with such an 
additional variable. 


7.3 Average Mach Number 

The eigenvalues of the one-dimensional average jacobian (44) are obtained from 
(105)-(106) as 

Ai = ? (112) 

P 

A 2 ,3 = -X + u d ± {u\ + a 2 ) 1 (113) 

P v ' 

Expression (113) leads to an average-Mach number definition as follows. For a 
positive velocity, the eigenvalues (112)-(113) are all non-negative if 


which yields 


— + u d > 
P 




+ 2 —u d > a 2 

P 


(114) 


(115) 
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Consequently, the average eigenvalues have all uniform sign if the expression 


M = 



(116) 


is greater than one. Therefore, M defines the average Mach number. The same 
result is achieved by enforcing non positivity of A 2 , for a negative velocity. Expres- 
sion (116) depends on both q and the partial derivatives of pressure, and in general 
differs from the classical expression evaluated at the ideal-gas average state, unless 
the average partial derivatives correspond to q itself, see Section 8. 


7.4 Average Eigenvectors 

The three right eigenvectors for the one-dimensional similarity transformation 


- (?£ 

9,p, U, 


= XAX 


-1 


(117) 


are 



(118) 


where the second and third column for X respectively correspond to the positive 
and negative sign of the ensemble ± in the second column of (118). Therefore, the 
average conservation-variable eigenvalues, Mach number, and eigenvectors (105)- 
(106), (116), and (118) differ from the classical forms due to their dependence 
upon the average deviation velocity «<*, in addition to the average sound speed 
squared a 2 . Consequently, the sign of (106), (113) may not coincide with the sign of 
the pointwise eigenvalues in transonic flow regions, while X may not revert to the 
pointwise eigenvectors. Hence, in equilibrium real gas situations, any approximation 
to these exact internal conservation variable average eigenvalues and eigenvectors 
may yield for some thermodynamic states a terminal flux splitting scheme with 
locally inadequate intrinsic dissipation. For an ideal-gas specification the average 
partial derivatives correspond to q itself, see Section 8. Hence, the deviation speed 
vanishes and consequently (105)-(106), (116), and (118) revert to the classical forms. 
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8 Partial Derivatives of Pressure 

In the published procedures [4,8,9] the average partial derivatives of pressure are es- 
sentially regarded as additional independent variables. Hence, they are algebraically 
determined using additional geometric projections and scale factors. In distinction, 
these derivatives are herein simultaneously computed by determining their intrinsic 
acoustic-wave average state, for the mean value theorem resolution of the pressure 
variation jacobian. 

8.1 Pressure Jacobian and Compatibility Relations 

For a homogeneous gas in thermo-chemical equilibrium, the pressure equation of 
state may be expressed as the function differentiable almost everywhere [l] in 

P = p(p, £ ) ( 119 ) 

where e denotes the mass-specific internal energy. Therefore, for a given thermody- 
namic state (p, e), both p and its partial derivatives are theoretically determined. 
For a CFD simulation, hence for a given < 7 , the mass-specific internal energy is cast 
as 

Therefore, s is no longer an independent variable, but intrinsically depends on 
the conservation state variable q. Consequently, within an equilibrium-gas CFD 
algorithm, p is expressed as 

P = p(p,^(p,m,E)) = p(q(x,t)) (121) 

Hence, the pressure jacobian is also a function of q 

Consequently, the average pressure jacobian corresponds to 

!■!<» <“> 

where q denotes the associated intrinsic evaluation state. 

The partial derivatives of pressure with respect to the state variable q are ex- 
pressed by way of the differentiation chain rule and the thermodynamic derivatives 
as 

(r) -(jr) +(jr) (r) =(¥) + (r) 

\ d P/m,E \ d PL \^ £ J p \ d Pj m,E \ d P I* \ d£ Jp P \ P J 

(124) 
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(125) 


dp \ _ ( dp\ / de \ _ ( dp\ ( m T \ 

d™/ P ,E~\ d£ )p\ d 7 n )r,E~\ d£ )p \ P 2 ) 



Whereupon these formulae lead to the fundamental compatibility relation 


(126) 



= 0 


pjn 


(127) 


The developments of Sections 4-5 yield the jacobian pressure variation expression 


R 

P ~P 



E l ) (128) 


Whereupon, insertion of the compatibility relation (127) into (128) and noting the 
dependence on the intrinsic evaluation state q yield the reduced form 


E r — E L ^j — -- • (rn R — 

(129) 

which is universally valid for both ideal and real gas equations of state, it only 
depends upon two partial derivatives, and it thereby closely resembles the classical 
thermodynamic expression 



f d£ 


p R -p L = m Ap+i^i a. 


'd P ' 

~de 


(130) 


Expression (129) coincides with a special form of (130). This result is obtained by 
inserting (124) and (126) into (129), which yields 


R 

P ~P 




'd P ' 
de , 


1 ( m T m rn AE 

— E 1 Ap — ■ Am + 


P 2 


This relation reverts to (130) provided that 


Ae(q,Aq) = 


' - — T ~ 

m ■ m 


p 2 


. m 
E j Ap — — • 


P 


AE 

Am d — — 
P‘ P 


(131) 


(132) 


This expression corresponds to the variation of (120) in terms of its first differential, 
and in general does not coincide with the simple independent variation 


Ae = e (§ H ) - e ( q L ) 


(133) 
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even though the difference between (132) and (133) is of order O (A q 2 ). When the 
acoustic-wave average state q coincides with the convection average state q\ then 
(132) coincides with (133) 


Ae(q,Aq) = 


f —T — 

m • m 


#) Ap- ™ - Am + ^ = e (? H ) - e (? L ) (134) 


P 


However, the convection average state q may not always satisfy (129) or (130). 
Therefore, within a CFD algorithm, the pressure variation (130) is in general ex- 
pressed as 

dp' 


P R ~P L = 


, d P, 


Ap + 



Ae ( q , A q) 


For a non-vanishing simple variation Ae, expression (135) is cast as 


(135) 






(136) 


which formally reverts to (130), and corresponds to the actual pressure variation 
relation utilized in several reported procedures, as detailed in Section 9. Following 
(122)-(123), relations (129) and (135) are properly expressed as 


P 


R 



P 


R 





Ac(q,Aq) 


(137) 

(138) 


Therefore, all of these derivatives can be simultaneously determined by identifying 
the single state q , as detailed in the next section. 


8.2 Intrinsic Evaluation State Determination 

The existence of an internal state q for the evaluation of the partial derivatives in 
(137) is ensured by the multi-dimensional mean value theorem (11), since p is a 
differentiable function with continuous partial derivatives almost everywhere in $7 
[1]. The determination of the intrinsic state q = q L + a (c[ R — q L ^j state is then 
achieved by the equivalent direct computation of the single weight coefficient a. To 
this end, (137) is re-expressed as 
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F( 5 )=g(/ + 5 (,''-/))) ra A, + 

(rn L + a ( m R - m L )) ^ 


AE 


P,m 


p L + a ( p R — p L ) 


• Am 


/ 


— A p = 0 

(139) 

This relation constitutes a single equation for the unknown 5 which can be solved 
to an arbitrary degree of accuracy by way of a Newton-Raphson iteration strategy 

(140) 


a s+1 — a’ — 


(—) 

\ Aa / a * 


where superscript “s” denotes the iteration index. Owing to the multi- dimensional 
mean value theorem, the solution 5 satisfies the inequality 


0.0 < 5 < 1.0 


(141) 


This relation circumscribes the root of (139), and an initial estimate a 0 for a first 
mesh-node computation may then coincide with (62), while Aa 1.0. Whereupon, 
a 0 at subsequent mesh nodes may coincide with the value of a at a neighboring point 
where (139) is already solved. This choice leads to an efficient computational strat- 
egy that only requires two or three iterations owing to the good initial estimate. 
Therefore, evaluating (140) is relatively inexpensive, it leads to a super linear con- 
vergent process, and directly yields a. Hence, the average state q as well as the 
average partial derivatives that satisfy (137)-(138) are simultaneously determined. 
This is the procedure utilized for the benchmark test in Section 10. 

A viable expression for A F in the denominator of (140) is 

(AF)„, s (Aj>,)„. • A p + (Ape )„. ■ (aE - ^j- ■ Am) + 

( Ap~- — Am] • AmAa (142) 

P / a 3 V P / a* 

where the last term corresponds to the differential of the coefficient of A ps with 
respect to a. Furthermore, (A p p ) aS and (A PE) a s are 

(A Pp ) al = p P (q L + (a* + Aa) (q R - q L )) - Pp (q L + a* (q R - q L )) 

(A p E ) a s = Pe(q L + (<** + Aa) (q R - q L )) ~ Pe(q L + a* (? R ~ ? L )) (143) 

Note that at covergence, the solution 5 no longer depends upon (142)-(143). 
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8.3 Ideal Gas Partial Derivatives 


For an ideal gas, (139) is identically satisfied for a coinciding with (62). This 
conclusion follows from the state equation 


P = (7“ 1) \ E- 

and associated average partial derivatives 

( d P 


T > 

m m 


2p 




1) 


— T — 

m * m 


2? 


m 

drn J P ,E 7 1 £ 


The insertion of (144)-(145) into (139) yields 


1 / m T ■ m 

A E A 

2 V P . 


m T ■ rn rn 1 

— A p- — 
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(144) 



p,m 


Am + A E 


hence 


— T — — T 

rri m n m 

— — A p + 2— — • Am — A 

P P 


f T \ 

m • tii 


= 0 


(7-1) 

(145) 


(146) 


(147) 


This expression is identically satisfied by the average density and momentum com- 
ponents (56)-(57), since (147), for a two-dimensional formulation, constitutes the 
sum of the first relation in (80) and the second relation in (84), which are them- 
selves identically satisfied by the convection average state. The same conclusion is 
achieved for a three-dimensional formulation. 


9 Comparisons with Other Formulations 


The pressure variation and sound speed expressions reported in [4,8,9] are ana- 
lyzed and their underlying hypotheses identified. These derivations algebraically 
determine the partial derivatives in (138) without reference to the intrinsic state q. 
The various solutions for each partial derivative set generally do not coincide with 
one another. Furthermore, it is a-priori unknown whether some evaluation state q a 
exists in the domain of the analytical partial derivatives such that 



(148) 


where subscript a denotes algebraic determination. And even if such states existed, 
they may correspond to external conservation- variable states. 
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9.1 Pressure and Sound Speed Expressions 

The pressure and sound speed expressions associated with the current formulations 
are listed in Tables 1 and 2 along with the the thermodynamic relation for the sound 



The reported formulae [4,8,9] use the thermodynamic derivatives of pressure and 
some combinations of flow variables in the respective sound speed expressions. Sev- 
eral derivations also introduce an extra average mass specific internal energy. The 
average derivatives in these constructions are synthesized using the assumptions 
that the general jacobian pressure variation and sound speed squared can be reli- 
ably satisfied using the following constraints 

1. ideal-gas average state evaluation of selected expressions within the pressure 
variation 
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2. selection of simple integration paths and associated auxiliary state 

3. geometric solution projection with arbitrary scale factors 


Significantly, the developed exact expressions (111)-(129), which are independent 
of these assumptions, lead to the reported formulae when these constraints are 
explicitly used. 

The equation of state used by Glaister [4] and Liou et al. [8] is 


P = P(P,£) 

Thus, the partial derivatives of pressure in (44) are cast as 


m T ■ m 
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(150) 

(151) 
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The insertion of (150)-(152) into (129) yields 



Ae 


Ae(q, A q) 


( mi m e\ rrv 


A E 


(153) 

Whereupon, the intrinsic assumption in the formulations of Glaister and Liou et al. 
is that (153) can be reliably satisfied by both evaluating the terms between brackets 
at the ideal-gas average state, and introducing the average mass specific internal 
energy 

With this assumption, the components of the deviation velocity Ud vanish, while 
(153) becomes 


(154) 


Ap = 


'dp' 

, d P, 

'dp' 

, d Pj 
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which coincides with (136) and the expression used in [4,8]. The average partial 
derivatives in (155) are then determined by projection, as highlighted in Section 
9.2. Note that even if these partial derivatives were exactly determined, using the 
procedure of Section 8.2 for example, there would be no a-priori assurance that the 
evaluation state (p, e) corresponds to an internal conservation-variable state. The 
expression for the speed of sound in terms of p and e is then obtained by inserting 
(150)-(152) into (111) as 



Ae 


Ae(q,Aq) 


({ 


H 


~ m 
£ z- 


( m m\ 1 


(156) 


With the stated assumptions, (156) becomes 


( dp\ ( dp\ \H — £ fh T ■ m 

W, + yai)], [“ ? W~. 


(157) 


which coincides with the corresponding expression in [4,8]. Considering that (154) 
is not evaluated at an internal conservation- variable state, but results from a direct 
average, expressions (155)-(157) do no revert to the standard ideal-gas expressions. 
Finally, the internal-state average jacobian (44) is also modified by restating the 
partial derivatives of pressure therein, with expressions (150)-(152) with coefficients 
evaluated at the ideal gas average. This operation then yields the matrix utilized in 
[4,8]. Therefore, it is not a-priori guaranteed that this matrix is evaluated at internal 
conservation-variable states, with unpredictable consequences on the eigenvector 
matrix, hence on the accuracy of the associated flux-difference-splitting results. 

Similar conclusions carry over to the developments of Vinokur et al. [9]. Their 
reference equation of state is 

- P = P(P,pe ) ( 15S ) 

Consequently, by using (120) the partial derivatives of pressure in (129) are cast as 


(dp\ /dp\ ( 9p \ (rn^_ni\ 
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The insertion of (159)-(161) into (129) yields 
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Whereupon, the intrinsic assumption in the formulation of Vinokur et al. is that 
(162) can be reliably satisfied by evaluating the terms between brackets at the ideal- 
gas average state. With this assumption, the components of the deviation velocity 
Ud vanish, while (162) becomes 


A p 



Ap + 

p* 



AE- A 


f T > 

m • m 


2 P 



A p + 

PC 



(163) 


which coincides with the corresponding expression used in [9]. The average partial 
derivatives in (163) are then determined by projection, as highlighted in Section 
9.2. The same conclusions stated after (155) also apply to (163). The expression 
for the speed of sound in terms of p and pe is then obtained by inserting (159)-(161) 
into (111) as 


c 


2 





(164) 


With the stated assumption, (164) becomes 
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2 




— T — 

m • m 

2 p 2 


(165) 


which coincides with the expression in [9]. Expressions (163)-(165) then revert to 
the standard ideal-gas expressions, since no additional average mass specific internal 
energy is utilized. Finally, the internal state average jacobian (44) is also modified 
by restating the partial derivatives of pressure therein, with expressions (159)-(161) 
with coefficients evaluated at the ideal gas average. This operation then yields 
the matrix utilized in [9], which again is not guaranteed to depend consistently on 
internal conservation-variable states. 


9*2 Algebraic Partial Derivatives 

The average partial derivatives in (155), (157), (163), and (165) are algebraically de- 
termined without specific reference to their intrinsic evaluation state. This algebraic 
determination accrues from a projection of a set of trial partial derivatives onto the 
admissible-state line. These trial derivatives are obtained through an approximate 
evaluation of the line integrals in the exact pressure variation expression 
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Therefore, the trial partial derivatives are defined as 

J_ f (dp 

A p Jr \dp 




dp , hf * 



(166) 


(167) 


while (166) is one form of the admissible-state line equation. Note that the devel- 
opments of Vinokur et al. are obtained from (166)-(167) by replacing e with pe. 
At this point, a multitude of integration paths can be selected for the exact inte- 
gration of dp since this is an exact differential and consequently A p is independent 
of the integration path. Nevertheless, the individual partial derivatives in (166) do 
depend upon the integration path, even though A p does not. Hence, different sets 
of partial derivatives (167) and associated speeds of sound are in general obtained 
with different paths, which in turn will depend upon selected auxiliary states. The 
specific auxiliary states employed in [4,8,9] are 

!=(;,!*) , B = (At") , c = (?,?) (168) 

while the selected integration paths are 


Table 3: Integration Paths 


Author 

r 

Glaister 

[4] 

LABR, LBAR 

Liou et ah 

[3] 

LCR 

Vinokur et al. 

[9] 

LR 


With these paths, the chosen forms of (167) are 


Table 4a-b: Trial Partial Derivatives 


Author 


f 

w £ 

Glaister 

Liou et al. 
Vinokur et al. 

Approximal 

P A ~P L 

Ap 

0“ 

)e integration of (167) 

Author 


l de ) 

A 

Glaister 
Liou et al. 
Vinokur et al. 

P^P L 

(#■0 

Approximate integration of (167) 
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The projections of these approximations onto the admissible- state straight line, 
corresponding to (166), are then the coordinates of the point shared by this line and 
the approximation line. The latter line contains the trial derivative approximation 
point, with coordinates based on (167), and intersects the admissible-state line at a 
suitable angle, as determined by selected scale factors. Representative closed form 
solutions are 


Table 5a-b: Average Partial Derivatives 


Author 


( ^p\ 

\&L 


Ap , p A - p a 

A p A p J 


Glaister [4] 
Liou et al. [8] 


A p 
A p 


_ (??\ — + (^e\ 

\ d£ ) p A P \ d P) e 


Author 





Glaister 
Liou et al. 


A p p~ - p' 
Ae + As 


A p _ fdp\ A p fdp\ 
Ae \ d Pj e Ae \de ) 


where Glaister’s trial derivatives from Table 4 are used to obtain his solution from 
Liou’s determination in this table. 

Analytically, these average partial derivatives correspond to the solution (x, y) 
of the linear system 


j kx + f3y = 7 
( (3x — ny — f}x — KtJ 


(169) 


where the first equation is associated with (166), while the second represents the 
approximation line. The solution of (169) is 


+ 0 {fix - «y) 
k 2 + (30 

07 - K ( 0 x - Ky'j 

V = K 2 + 00 

By author, the variables in (169)-(170) are defined as 


(170) 
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Table 6a-c: Variables in (169)-(170). 


Author 

X 

y 

Glaister [4] 

KVA/. 

l ( 

£\ d P/c 

Liou et al. [8] 

1 (dp\ 

1 / dp\ 

£\ d P/s 

Vinokur et al. [9] 

(I)." 

( dp\ (dp \ 
\dp) m \dpe) 0 


Author 

X 

y 

Glaister 

1 ( dp'" 
7l\de, 

\ 1 ( 

L 

Liou et al. 

1 (dp' 
71 \c?e, 

\ 1 ( dp \ 

K £ \ 9 p). 

Vinokur et al. 

uy, 

(dp\ ( dp\ 

, \Sp)„\ 9 Pt), 


Author 

K 0 

0 

7 

Glaister 

_ Ae _A p 
71— S-r- 


1 

Liou et al. 

A p A p 

_ Ae . Ap 
71— £-r~ 

Ap 

£ ir 

1 

Vinokur et al. 

A p A p 

A p —Ap 

Ap 

-S 2 Ap 

A (pe) 


Note that since Glaister and Liou et al. select 0 = 0, the admissible-state line is 
orthogonal to the approximation line in their implementation, whereas this feature 
is not present in the formulation of Vinokur and Montagne. The parameters /c, 
0, and 0 depend upon 7 Z, £, and S which are free are non-dimensionalizing scale 
factors. Therefore, even after specifying an integration path, a formulational non- 
uniqueness still persists upon variation of these scale factors. 

With the specific selection 




(171) 


and using the trial derivatives in Table 4, expressions (170) yield the solutions of 
Glaister and Liou et al. listed in Table 5. All of these average partial derivatives 
algebraically satisfy the pressure relation (166), even though they may become 
indeterminate for vanishing A p and Ae, and generally do not coincide with one 
another. 


37 


















10 Representative Results 

As indicated in Section 2, the solution of (14) and (139) is in general non-unique. 
This feature is now demonstrated for a shock-tube test case [4] for which intrinsi- 
cally exact average states are determined without separating the kinematic/energy 
terms from the pressure relations. With these solutions, the associated exact ther- 
modynamic pressure derivatives and average sound speed squared are determined 
as indicated in Section 8.2. These results are then compared with the predictions 
of Liou’s and Glaister’s methods for an appropriate quantification of their inher- 
ent approximations, as induced by the algebraic synthesis of the average partial 
derivatives of pressure. 

The test left and right states are 


Table 7: 


Initial Shock Tube S tates 



Left 

Right 

p 

100 atm 

1 atm 

T 

9014 K 

300 K 

u 

0 m/sec 

0 m/sec 


Using the left-state pressure and temperature as reference values, the associated 
non-dimensional states are 


Table 8: Non-Dimensional States 



Left 

Right 

A 

B 

P 

1.000 

0.447 

0.447 

1.000 

m 

0.000 

0.000 

0.000 

0.000 

E 

5.491 

2.498E-02 

2.454 

5.580E-02 

£ 

5.491 

5.580E-02 

5.491 

5.580E-02 

p 

1.000 

1.000E-02 

0.403 

2.231E-02 

T 

1.000 

3.327E-02 

0.948 

3.326E-02 


1.053 

2.231E-02 

1.008 

2.2311E-02 

\ d p)e 





/ op \ 

0.140 

0.178 

5.880E-02 

0.399 

\de ) 

K ' p 
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Figure 1 shows the variation of (139) for 0.0 < a < 1.0, for both ideal and real 
neutral air, corresponding to the left and right states in Table 8. 



Figure 1: Variation of Function (139) 

Evidently, the ideal-gas version of (139) is identically satisfied for any value of a, 
following (146)-(147), since m = 0. Conversely, for real air, not only does a solution 
exist, but as many as three are possible in this case. These solutions, as generated 
by (140) with Aa = 1.0 E — 04, are listed in the following table, along with the 
ideal-gas average state 
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Table 9: Ideal-Gas Average State and Exact Solutions. 


a 

P 

m 

E 

e 

Ae(q,Aq) 

Ae 

Aejq, A q) 

As 



Average 

Sol. I 

Sol. II 

Sol. Ill 

0.599 

0.689 

0.860 

0.908 

0.669 

0.618 

0.524 

0.498 

0.000 

0.000 

0.000 

0.000 

2.216 

1.719 

0.789 

0.527 

3.312 

2.782 

1.505 

1.058 

-5.435 

-6.356 

-8.842 

-9.800 

-5.435 

-5.435 

-5.435 

-5.435 

1.000 

1.169 

1.627 

1.803 

0.709 

0.618 

0.369 

0.290 

9.712E-02 

0.102 

8.895E-02 

8.310E-02 

9.712E-02 

0.119 

0.145 

0.150 

0.797 

0.806 

0.778 

0.767 


According to these results, the variation of mass specific internal energy A e(q, A q) 
significantly depends upon the intrinsic state q. Conversely, the dependence of the 
sound speed squared on q is less pronounced, while the eigenvalues corresponding 
to these three solutions all display the same algebraic sign structure, in agreement 
with the eigenvalue consistency theorem. 

For this test case, the trial and projected average pressure derivatives for Glais- 
ter’s and Liou’s methods are 


Table 10: Trial and Average Partial Derivatives 





^ 


1 

Author 


/ 

\ d PJ * 

( 2 ). 

( 2 ), 

( 2 ), 

Glaister 

[4] 

■ESI 

0.551 

0.180 

0.126 

Liou et al. 

[8] 


0.772 

9.712E-02 

0.103 
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Furthermore, the admissible-state straight-line equation associated with (166) 
and the first equation in (169) is 

(|), = - 0102 ( l )/ 0 - 182 (172) 

Figure 2 graphs the admissible- state straight line along with Glaister’s and Liou’s 
solution projections. 



Figure 2: Admissible-State Straight Line 

The trial solution of Glaister’s method is farther than Liou’s from the admissible 
state line, even though its associated projected solution is closer to all three exact 
solutions than Liou’s. Clearly, a different choice of the scale factors 71 and S may 
yield closer terminal projections, even though such factors would remain solution 
dependent, and there seems to exist no practical way to determine them accurately. 

The solution closest to the ideal-gas average state is then isolated, since with 
this selection the exact average eigenvalues, Mach number, and eigenvectors (105), 
(106), (116), and (118) converge to the pointwise forms. Table 11 presents the exact 
solution I and the ideal-gas averages, as compared with the predictions of selected 
expressions in Tables 2 and 5. 
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Tab] 

[ell: Exact and Approximate Solutions 


Exact 

Average 

Glaister 

Liou et. al 

( ^P\ 

W. 

e% 

0.618 
0.00 % 

0.709 

-14.56% 

0.551 
-10.84 % 

0.772 
+24.91 % 

( 

[Ts 

0.119 

9.712E-02 

0.126 

0.103 

\ / p 

e % 

0.00 % 

18.38% 

-5.88 % 

+13.44 % 

c 2 

0.806 

0.797 

0.666 

0.867 

e% 

0.00 % 

1.11% 

-17.36 % 

+7.56% 


For the thermodynamic partial derivatives of pressure, the percent error e% 
reaches 24.91%, while that incurred in the speed of sound squared is as high as 
17.36 %. Furthermore, in this case, these errors are not significantly different than 
those generated by simply evaluating all thermodynamic partial derivatives at the 
ideal-gas average state, as indicated in Table 11. In general, these errors will prop- 
agate to the eigenvector matrix with unpredictable consequences on the accurate 
flux- difference splitting computations of subtle flow features such as shock waves 
and contact discontinuities. 

Concluding Remarks 

The theoretical developments in this report yield unambiguous internal-average 
conservation variables for curvilinear- coordinate multi- dimensional flux difference 
resolution, including the speed of sound and the partial derivatives of pressure. 
These derivatives are determined for arbitrary equations of state without any extra 
assumption. Hence, no average internal energy is introduced, and no approxima- 
tions and projections are employed. 

The Jacobian matrices (44) and (71) attain a tidy structure, and for equilib- 
rium real-gas equations of state, their eigenvalues and associated Mach number are 
shown to differ from the classical forms. Only after introducing subtle assumptions, 
documented in Section 9, can the average eigenvalues formally coincide in these con- 
ditions with the classical forms. Of course, the unavoidable theoretical price to pay 
is the loss of an a-priori assurance that the Jacobian partial derivatives of pressure 
and speed of sound either exist in the domain of the analytical partial derivatives, 
or correspond to an evaluation at internal-average conservation- variable states, with 
unpredictable consequences on the eigenvector matrix, and hence on highly accurate 
flux-difference-splitting computations of subtle flow features like shock waves and 
contact ^discontinuities. Conversely, the formulae presented in this research afford 
this a-priori assurance. Furthermore, their very existence, uniqueness, and numer- 
ical demonstration lead to the fundamental conclusion that the generalization of 
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the average-state notion to an equilibrium real-gas equation of state introduces no 
additional degree of freedom in the determination of the average partial derivatives 
of pressure, unlike what is stated in [9]. 

Certainly, the reported average jacobians will converge to their respective classi- 
cal pointwise forms when the given computational-cell left and right states approach 
each other, which for smooth flows occurs under appropriate mesh refinements. 
Therefore, although these average jacobians may correspond to external average 
conservation variables, with thermodynamic partial derivatives and speed of sound 
possibly affected by non negligible errors at shock locations, their utilization in sig- 
nificantly refined grid computations, e.g. [8], has led to satisfactory results. One 
however wonders whether these procedures depend upon uniform refined grids for 
their reliable operation. 

The developed average jacobians and associated eigenvalues and eigenvectors, 
consistently depend on internal- average conservation variables and without extra 
assumptions they exactly satisfy the flux difference resolution even for arbitrary 
separation of the computational-cell left and right states. Furthermore, they are 
valid for arbitrary multi-dimensional curvilinear-coordinate meshes. 
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